clc;clear;
tic;
format short
% Pauli matrices
ai=sqrt(-1);

sx=[0 0.5; 0.5 0];                  % matrix of S_x for spin 1/2
sy=[0 -0.5*ai; 0.5*ai 0];           % matrix of S_y for spin 1/2
sz=[0.5 0; 0 -0.5];                 % matrix of S_z for spin 1/2
od=[1 0; 0 1];                      % 2*2 unity matrix

unit=eye(4);

sz1=1;sz2=1;
sz1=kron(sz,od);
sz2=kron(od,sz);

sx1=kron(sx,od);
sx2=kron(od,sx);

sy1=kron(sy,od);
sy2=kron(od,sy);

rho_ortho=3/4*unit+(sx1*sx2+sy1*sy2+sz1*sz2)
rho_para=1/4*unit-(sx1*sx2+sy1*sy2+sz1*sz2)
